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ABSTRACT 

We study the accumulation of errors in cosmological A^-body algorithms that are 
caused by representing the continuous distribution of matter by massive particles, 
comparing the P'^M and Adaptive Multigrid codes. We use for this a new measure 
of two-body relaxation suitable for nonstationary gravitating systems. This measure 
is based on accumulated deflection angles of particle orbits and it does not saturate 
during evolution. We have found that the role of gravitational collisions in standard 
P^M-models is rather high, and in order to avoid collision effects we recommend to 
use comoving softening parameters e >— 1.0 (in grid units). 
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1 INTRODUCTION 

Cosmological simulations use A'^-body algorithms that re- 
place the initially continuous distribution of (dark) matter 
by discrete and rather massive 'particles'. This helps to sim- 
ulate intersecting streams of coUisionless matter that are 
common at the time of formation of cosmological structure. 
Interactions between these particles may, however, differ 
from the dynamics of the continuous media they are meant 
to approximate. The main source of discrepancy are close 
encounters between the particles (gravitational collisions). 

The popular belief has been that if the number of points 
in the simulation is sufficiently large, the discreteness effects 
(gravitational collisions or close encounters between mass 
points) do not play an important role. These effects have 
been studied only for tree-codes in case of stationary grav- 
itating systems (Hernquist & Barnes 1990, Huang, Dubin- 
ski & Carlberg 1993). We shall describe in this paper the 
role of gravitational collisions in high-resolution cosmologi- 
cal codes. 

Among the different N-body algorithms used for cos- 
mological simulations the P'^M (Particle-Particle-Particle- 
Mesh) algorithm has become the industry standard for high- 
resolution simulations. It was designed initially for plasma 
physics simulations and ported to cosmology by Efstathiou 
& Eastwood (1981). It has been clearly documented in text- 
books (Hockney & Eastwood, 1988, Couchman, 1995) and 
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the code is practically in public domain. This code was the 
first candidate for our study. 

The P^M-code improves considerably the resolution of 
a standard PM (Particle-Mesh) code by splitting forces in 
two parts - the long-range force from the overall density dis- 
tribution on a mesh and the short-range force from nearest 
neighbours, calculating the latter by summation over pairs 
of particles. The algorithm was devised initially to simu- 
late collections of pointlike particles. In cosmological simula- 
tions of formation of structure the basic dark matter density 
is smooth and does not lend itself easily to representation 
by discrete mass points. This difficulty is usually alleviated 
by introducing smoothing parameters in pairwise forces, ef- 
fectively replacing the mass points by extended spherical 
clouds. 

In order to study dicreteness properties of the P'^M- 
code one should have smooth codes with spatial resolution 
comparable to the P'^M. Such codes have appeared only 
recently, and most of them are grid-based multi-resolution 
schemes (except that of Pen (1994) that uses a global de- 
formed grid). These codes improve the spatial resolution in 
selected regions (mostly in those of high density). We shall 
describe these codes and their main differences below; for 
comparison we used our adaptive multigrid code (Suisalu & 
Saar 1995). 

We ran also one simulation using the popular tree- 
code (Barnes & Hut 1986). Its deflection properties have 
been studied before (Hernquist & Barnes 1990, Huang et al. 
1993), but for a stationary case only. 

The usual tool to measure pairwise collisions is the 
study of energy diffusion (see the studies of three-code sim- 
ulations cited above) . This is proper for stationary gravitat- 
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ing systems, where energy is otherwise conserved and heat- 
ing comes only from pairwise interactions; in cosmology the 
overall energy changes with time and it is not the best mea- 
sure for collisions. We propose instead of it to study orbital 
deflection angles, concentrating only on the change of direc- 
tion of velocities and not on the change of their absolute 
values. 



2 SIMULATIONS 

As the adaptive grid methods are not so well known as the 
P^M and one is forced to select between different versions, 
we shall describe their differences below. 

AU these methods form subgrids of higher spatial reso- 
lution in regions where density is higher than a prescribed 
limit. They use for boundary conditions on a subgrid the val- 
ues of the potential (Dirichlet conditions) interpolated from 
the coarser grid, but they differ when finding the solutions 
for the potential. The problem is if the coarse grid solution 
should depend on the finer grid solution at the same point, 
or not. It certainly should if we used exact direct solvers 
of linear systems arising from discretized partial differential 
equations, but the situation is not clear for approximate it- 
erative solvers. In cosmological simulations we do not have 
an exact source term, either, as our density determination 
may depend on a grid level. 

In an ideal situation solutions of different resolution 
should converge rapidly during the iteration process. In 
practice convergence is rather slow, as solutions on coarser 
grids are considerably changed by temporary solutions on 
finer subgrids (see, e.g. Brandt 1987). This issue has been 
discussed by Anninos, Norman & Clarke (1994) who clas- 
sify adaptive grid methods as having either one or two-way 
interfaces between parent and subgrids. 

As an example, in the multi-grid method developed by 
Jessop, Duncan & Smith (1994) only one-way grid interfaces 
are used that implies that the local fine grid solution can be 
different from the coarse grid solution. They find the local 
solution by iteration on the subgrid only. A recent adaptive 
code by Splinter (1995) uses a similar methodology, although 
the algorithms for finding the solution on subgrids difi'er. 

In contrast, in our AMG (Adaptive Multigrid) code (Su- 
isalu & Saar 1995) the solution for local subgrids is obtained 
simultaneously with the global solution using the full multi- 
grid method (Brandt 1977). According to the above clas- 
sification AMG uses two-way grid interfaces between finer 
and coarser grids, as information passes in both directions 
during iterations. The reason why we have been left alone 
in this class is that here it is harder to get good convergence 
of the iterative solution process. We believe that two-way 
interfaces, once they have been built, are closer to the exact 
solution, and we shall use below our AMG-code for compar- 
ison with the AP^M-code. 

At first we encountered difficulties in building our mod- 
els, as the natural boundary conditions for a multigrid code 
are the Dirichlet conditions, and our code was initially tuned 
for this case. In order to follow the evolution of spatially pe- 
riodic P^M models we had to modify our code for periodic 
boundary conditions. 

If we use iterative methods to solve the Poisson equation 
on a grid, we discover that in the case of periodic boundaries 



Table 1. Simulation parameters 



Sim. 


L 


N 


e 




P3M1* 


32^ 


323 


0.2 


1 


P3M2 


32^ 


323 


1.0 


1 


P3M3 


32^ 


323 


2.0 


1 


P3M4 


32^ 


643 


0.2 


1 


APM 


32^ 


323 


0.21' 


1/8 


AMGl 


32^ 


323 




1/16 


AMG2 


32^ 


643 




1/32 


TREE 




16440 


0.2 





* Three different spectra 

t Short-range accelerations excluded 



the linearized system of algebraic equations that we have to 
solve becomes singular. Singular systems are more difficult 
to solve, and the problem becomes even harder if we con- 
sider interactions between the subgrids and the global grid. 
In the case of full two-way interfaces that we use, solutions 
on subgrids always induce changes in global solutions that 
violate global boundary conditions, and this requires usually 
additional iterations. In the periodic case the global solution 
is very sensitive to violations of the zero total mass condi- 
tion that are generated by local subgrids, and convergence 
becomes very slow. 

The remedy we chose was to use the Galerkin condition 
to modify the difference operators on coarser grids. Namely, 
we calculate the difference operator on the coarse grid 
as 

L" = 4^L'^/^, (1) 

where L*^ is the difference operator on the fine grid, is 
the restriction operator from the fine grid to the coarse grid 
and is the corresponding interpolation operator. This 
technique helps to integrate singular systems, and it is de- 
scribed in detail by Hackbusch (1985). We found it essential 
for speeding up the convergence of the solution for the grav- 
itational potential. 

We chose to impose these conditions only when find- 
ing the solution for the global grid and did not modify the 
differential operator L on subgrids. In order to satisfy the 
boundary conditions themselves, we copied appropriately 
the boundary regions on global grids between iterations. 

For p3M simulations we used H. Couchman's Adaptive 
p3M code Ap3M (Couchman 1991) that speeds up the nor- 
mal p3M by generating subgrids in regions of high particle 
density and finding there a smooth solution for the poten- 
tial. This works to decrease the volume of pairwise force 
summation and considerably speeds up the algorithm. As 
concerns the adaptive smooth solution, the Ap3M-method 
belongs to the class of those with an one-way interface. 

We have run and analyzed 6 p3M simulations, 3 PM- 
type simulations and one tree-code simulation, their param- 
eters are summarized in Table 0. The first column labels the 
simulations, L is the base mesh size in cells (the physical 
size is SOh^^ Mpc for all simulations), A'^ is the number of 
particles used, e is the comoving softening parameter and 
hmin the minimal meshsize for adaptive models. 

Couchman (1995) advises using softening radii that are 
constant in physical coordinates that leads to increasing dis- 
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Figure 2. Particle distribution for the model AMGl (AMG with 
4 subgrid levels) at a = 9. The initial conditions and the cube size 
are the same as for the model P3M1 in Fig. |^. 

creteness during simulations. While this could be appropri- 
ate to decribe the interaction of separate mass concentra- 
tions that may form during the simulation, it is certainly 
not the best description for dark matter. We used comov- 
ing softening to keep P^M-simulations closer to smooth grid 
simulations; e.g., Gelb & Bertschinger (1994) have chosen 
the same approach. 



Altogether we used three different values of the soft- 
ening parameter e = 0.2, 1 and 2 (our minimal softening 
is that normally used in P"^M-simulations). The maximum 
number of refinement levels for AP"^M was 4. The model we 
call APM is similar to the P3M1 with the only difference 
that we commented out in the source code the lines that 
updated velocities by accelerations from short-range forces 
- this is essentially a PM-model with adaptive mesh refine- 
ment, hidden inside the camp of AP"^M-models. 

Two other PM simulations labeled AMGl and AMG2 
were run using our Adaptive Multigrid code modified for 
periodic boundary conditions as described above. This code 
gives a spatial resolution for forces similar to that of P^M 
codes without the need to consider pairwise forces. The last 
two codes differ by the number of points used and also by 
the number of subgrid levels allowed, 4 for AMGl and 5 for 
AMG2. 

The tree-code we used was the so-called 'Barnes' export' 
code with quadrupole corrections. While the computational 
volume for all other models was a periodic cube with the 
size of 80/i~^ Mpc, the volume for the tree-code was a sphere 
(with the same diameter). We used an opening angle 6 = 0.8 
and physical coordinates with the timestep of 0.001 Gy. 

We used the same initial conditions for all our models, 
trying to eliminate all possible sources of differences. The 
initial conditions were generated using the test power law 
spectrum with the spectral index n = — 1 from the Couch- 
man's Adaptive P^M distribution (see Couchman 1995 for 
a detailed description). For the P3Ml-model we generated 
also two other realizations with the initial spectra as power 
laws with indices n = —2 and n = 0. When we compare 
different models from the Table 0, we use always those with 
the n = — 1 spectrum. 

As we simulated the simple Q, = 1 cosmology, we could 
choose our starting time at will (we use the scalefactor a 
as our time variable). We started at a = 1 and followed 
the evolution until a = 9 when erg ~ 1 for all models, thus 
we brought our models up to the present time. In order to 
study the growth of orbital deflections further, we followed 
the evolution of selected models for longer periods, up to 
a = 30. 

We did not use the usual energy condition to check our 
time steps. Instead of this we use in our AMG code the 
Courant condition: 

(2) 

(we change the time step when necessary, using an asymmet- 
ric leapfrog integration). As the AP'^M-code uses a constant 
timestep, we ran first our AMGl-model, found the minimum 
timestep used (Aa = 0.0625) and used it for all simulations 
(we even ran AMGl once more). This guarantees the use of 
the same algorithm for integration in time and eliminates 
one possible source of difference between models. The value 
of Aa we chose corresponds to Ap — 0.01 for the time vari- 
able p = (3/2)a used in AP^M for Q ^ 1. This time step is 
much smaller than usual, but it is necessary to follow accu- 
rately particle trajectories. 

The typical density distribution for our models is shown 
in Figs. |l| and ^. The first one shows the particle distribu- 
tion for the model P3M1 at a = 9 and the second one a 
similar one for the model AMGl. The density distribution 
in the former seems to be more developed and has more dis- 
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Figure 3. The evolution of the initial density power spectrum 
P{k) (the curve marked a = 1) in different models. The spectra 
for the three models shown at the final moment a = 9 are labeled 
according to Table |l[ The wavenumber k is given in the comoving 
grid units. 



tinct substructure, while AMG has retained more of linear 
structure elements and is less clumpy. We may suspect that 
the difference is due to pairwise collisions but we cannot say 
this on the basis of comparing density distributions only. 
The overall impression of both models is rather similar, of 
course - they start from identical initial conditions. 

The evolution of the spectrum for different models was 
also rather similar (see Fig. The difference is in the 
growth of small-scale modes that is the highest for the P'^M- 
models with a small softening (P3M1), the lowest for the 
APM model (this is essentially a pure PM-model) and inter- 
mediate for the smooth but high-resolution AMGl run. The 
difference between the last two models can be explained by 
different depths of adaptive subgrid levels used (3 for APM 
and 4 for AMGl) and by the fact that the effective smooth- 
ing in P^M is somewhat higher than that in AMG (about 
3.5 versus 2 cell sizes). This makes the resolution achieved 
by AMGl about 4 times higher than that in APM. 



3 DEFLECTION ANGLES 

It is well known that during the linear stage of growth of 
density perturbations particle velocities retain their initial 
directions; the well-known Zeldovich approximation says 



X = q + b{t)v{q), 



(3) 



where x is the comoving coordinate of a particle labeled 
by its initial coordinates q and b{t) is the density growth 
rate from the linear theory, a function of time only. This 
shows clearly that while the velocity amplitude may change 
in time, its direction remains fixed. In fact, the components 
of the velocity direction are adiabatic invariants of motion 



in the initial stage, while particle energies change with time 
even without any interaction. This suggests that the change 
of the velocity direction is a better measure of gravitational 
interaction than the usually used energy diffusion. 

Nonlinear effects - growth of structure and pairwise col- 
lisions - both contribute to the change of the direction of ve- 
locity, and if we define the defiection angle dcj) for a particle 
by 



' = \v ■ dv\/v 



(4) 



then the total angle (p accumulated during the run of the 
simulation refiects all gravitational collisions suffered by the 
particle. The quantity (j) does not describe the real angle 
between the initial and final direction of the velocity of a 
particle, as it omits the other degree of freedom 6 necessary 
to describe a direction in a 3-D space. However, in case of a 
small total deflection the two angles are connected by 



(sin^ $) = 2/3[l - exp(-3S'/2)] 



(5) 



(Standish & Aksnes 1969), where the accumulated square 
deflection for a trajectory S is 



S = ^(d,^)^ 



(6) 



(a is our time coordinate). The mean in (^ concerns all 
trajectories with the same di^-sequence but with different 9. 
For small values of S it is equal to the mean square of the 
final deflection angle <E>, but for a long history of collisions 
(sin^ $) will reach the value 2/3 that describes an isotropic 
distribution. On the contrary, the accumulated square de- 
flection S does not saturate and continues to describe the 
total effect of gravitational interactions. In this respect it is 
also better than the often-used measure of orbital divergence 
that saturates easily for spatially limited systems. 

In order to be able to catch large single deflections we 
computed the defiection angles in the simulations as 



= arccos(t;t ■ Vt+st/\vt\\vt+st\), 



(7) 



this did not add significantly to the cost of the simulations. 
We summed both the deflection angles and their squares for 
every trajectory during a simulation. 

After this work was finished, we learned about a similar 
measure proposed recently by Bagla & Padmanabhan (1995) 
to describe nonlinearity in cosmological structure formation. 
They describe deviations from the linear stage by a measure 



Dgu = {u ~ gf ju 



(8) 



(as u is the velocity with respect to the dimensionless time 
a, its dimension is the same as that of the acceleration g). 
Their measure is easier to calculate than ours, but it has to 
be modified for other cosmologies (the Zeldovich approxima- 
tion gives zero for this expression only for f2 = 1 universe), 
ours does not depend on the background cosmology. But 
we see clearly that there are two sources for orbital defiec- 
tions - nonlinear evolution of particle orbits in the smooth 
background field described by ^ and the pairwise (gravi- 
tational) collisions we are looking for. Of course, both these 
sources contribute to the energy diffusion, too. 

While deflection angles reflect histories of individual 
particles, the simplest characteristic describing the change 
of all particle trajectories in a simulation is the mean of 
accumulated squares of deflection angles for all trajectories 
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Figure 4. Growth of the mean accumulated square deflection 
)S{ with dimensionless time a (the scale factor). The models are 
labeled according to Table |lj the highest growth is observed in 
case of the P^M-models with a normal smoothing parameter e = 
0.2. 



P3M1 




a 



Figure 5. Growth of the mean squared angular velocity de- 
flections )W{ with dimensionless time a. The models are labeled 
according to Table |^; the highest growth is observed in case of 
the P^M-models with a normal smoothing parameter e = 0.2. 



{S). If particle motions are quasi-linear, as in the Zeldovich 
approximation, particles follow their initial direction (only 
their velocity may change in time) and (5) = 0. Growth of 
the deflection angles describes nonlinear interactions, either 
via the mean fleld or by pairwise gravitational collisions. 

In Fig. ^ we see the evolution of the deflection measure 
{S) in time for different models. First we see a striking dif- 
ference between the P^M-models with different smoothing 
parameters; for the normally accepted smoothing parame- 
ter (e = 0.2) the accumulation of deflection angles is very 
rapid compared to those for larger e. Contrary to the com- 
mon belief, the rate of growth of deflection angles does not 
depend much of the total number of particles A'^ - compare 
the curves P3M1 and P3M4, the latter is for a model with 
8 times more particles than the P3M1. Our multigrid mod- 
els lie in the middle of the range, the model AMGl prac- 
tically coinciding with the P3M2. This is understandable, 
as e = 1 corresponds to a cell-sized smoothing; the model 
P3M3 e = 2 is evidently 'oversmoothed'. 

Surprisingly, the model AMG2 gives larger deflections 
than the AMGl. This can hardly be caused by the larger 
number of particles, a much more likely cause is that we 
have allowed an extra level of refinement in this model and 
follow particle trajectories better. An additional source of 
deflections could be the force errors that each refinement 
introduces near subgrid boundaries; as we have shown earlier 
(Suisalu and Saar 1995), these errors are usually small, less 
than one per cent, but in rare cases they may reach a few per 
cent. Similar errors are present in the AP'^M-code, where 
the usual requirement is to limit them by 6 per cent. We 
stress once more that this value describes very rare errors. 
Anyway, such errors may get amplified in a square- weighted 



characteristic as (S) is. Another surprise is a relatively low 
defiection measure obtained for our tree-code run. 

As we mentioned above, the deflection may be caused 
both by gravitational collisions and by interactions with the 
mean gravitational field. It is well known than in case of 
uncorrelated defiections particles follow a random walk in 
the energy space (see, e.g., Huang et al. 1993). This would 
translate in our case to 

^=n(A0)^ (9) 
da 

that would give a (5) ~ t dependence (n is here the mean 
frequency of collisions and {^<j)f the characteristic value of a 
single squared defiection). This is close to the a-dependence 
seen in Fig. ^. 

However, in our nonstationary simulations we have to 
consider also the possible role of the mean-field effects. We 
can roughly model strong mean-field deflections, supposing 
that all particles rotate in circular orbits with angular veloc- 
ities Ldi. In this case the computed deflection measure would 
grow as 

N 

i a i 

As we see, (S) caused by a strong mean-field interaction is 
also proportional to time a, as the expected effect of colli- 
sions was. 

In order to differentiate between the signatures of mean- 
field defiections and those caused by pairwise collisions, we 
introduced a new ('velocity') defiection measure W for the 
growth of the squared defiections of the angular velocity of 
a particle: 
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Figure 6. Growth of the velocity deflection measure {W) in the 
model P3M1 for different initial spectra with exponents n = 0, 
n = — 1 and n = —2. The dotted line shows a power law {W) ~ a. 



W{a) = 




where oj is the time-averaged rate of change of the deflection 
angle uj for a particle (the average angular velocity for planar 
orbits) during the whole simulation (from ao to a): 

Lu = l/(a-ao) / —da (12) 
/ da 

The velocity deflection measure W is, of course, zero for 
linear motions without the change of direction, but also for 
the case of strong interactions, when particles are trapped 
into systems and rotate with constant angular velocities uoi, 
<j)i = uJiU. In order to describe all particles we use a mean 
value of this measure (W). The growth of this measure in 
time for our simulations is shown in Fig. 

The difference of the mean velocity deflection measure 
between various models (Fig. ^ is considerably smaller than 
it was for the mean squared deflections (Fig. W). The P3M1- 
model leads the pack as before, with its many-particle ver- 
sion close behind, and the tree-code model has generated 
similar deflections to the high-resolution AMGl. All other 
models are relatively quiet. 

We studied the e = 0.2 case in more detail, trying to 
better understand the source of deflections. We continued 
the run more deep into the nonlinear times (Fig. |^. We 
can see that the growth, although rather rapid at the start 
of the simulation (this rise could be probably explained by 
non-self-consistent starting conditions), levels off to a power 
law with the exponent a « 1. As we have largely eliminated 
the mean-field effects, we may hope that this exponent tells 
us that the observed growth of {W) is caused by two-body 
effects. 

In order to check if this exponent depends on the large- 
scale environment we repeated the simulations for the model 




log W 

Figure 7. Distribution of the velocity deflection measure W for 
all particles in the P^M-model P3M1 for different moments a. 
The signature of high-deflection gravitational collisions is clearly 
seen as the right maximum. Observe how the distribution that 
is dominated by small deflections at the start of the simulation 
shifts gradually over to relatively large deflections. 



P3M1 for two more power spectra with the exponents n = 
—2 and n = 0. As can be seen from Fig. ^ the deflections 
grow in a similar fashion, and the exponents are practically 
the same - the curves differ by a multiplicative factor only. 
Examination of density distributions conflrms that structure 
is much stronger in the n — —2 case, the n = spectrum 
giving rise to a large number of small clumps and a rather 
diffusive sea of particles in between. This will naturally lead 
to a larger collision rate for the former model. 

It is clear that different particles follow extremely dif- 
ferent orbits, and in order to understand their evolution it 
is better to 'differentiate' the average quantitities shown in 
Figs. and to study the distributions of the velocity de- 
flection measure for our collection of trajectories. 

Fig. ^ demonstrates the evolution of this distribution 
in time for the coUisional model P3M1. The main feature 
of this distribution is the presence of two strong maxima 
at all times. We may suppose that the left maximum at 
smaller deflections is describing mean-fleld effects and the 
right maximum - gravitational collisions. During the evo- 
lution the number of particles that participate in collisions 
grows steadily. We see also the gradual shift of the distribu- 
tions towards larger deflections - the accumulated velocity 
deflections grow. The width of the right maximum that can 
be thought of as describing the number of strong collisions 
(about 3 per unit logarithmic interval) grows also. 

Similar distributions for the multigrid model AMGl are 
shown in Fig. ^ They are distinctly different from the P^M 
model shown before, having only one small-deflection maxi- 
mum. During evolution the distribution spreads and only the 
last distribution that corresponds to highly nonlinear stages 
of evolution of structure (more than two present lifetimes of 
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Figure 8. Distribution of the velocity deflection measure W for 
the multigrid model AMGl for different moments. The deflection 
measure grows gradually, but small deflections dominate at all 
times. 



Figure 9. Distribution of the velocity deflection measure W 
for the APM-model (without local forces) at different times. The 
distribution is initially smoother than for the multigrid models 
in Fig. ^, and the final highly nonlinear stage develops similar 
fluctuations. 



the Universe into the future) shows presence of an apprecia- 
ble high-deflection tail. This tail is caused by the growth of 
small-scale substructure by far beyond the present epoch. 

Even more smooth distributions can be found when run- 
ning the APM model, where the local pairwise forces were 
ignored. The distributions (Fig. |^ show only the signature 
of mean-field forces and do not evolve practically at all. This 
shows, first, a total absence of gravitational collisions typical 
for its parent P'^M-model, but also less substructure than 
in the AMGl-model (the spatial resolution of the present 
model was about 4 times lower than that of the AMGl). 
An exception is the distribution for large times that shows 
features similar to the AMGl, even with slightly larger am- 
plitudes. 

The tree-code shows features in between of the P3M1- 
model and the smooth models - the distribution of the ve- 
locity deflection measure W (Fig. [lo[ ) is rather wide, compa- 
rable to that of the model P3M1, but it has only one maxi- 
mum at all times. As the transformation between the phys- 
ical time coordinate t used here and the 'scale factor time' 
a used for other models is nonlinear, it is difficult to com- 
pare these distributions directly with the ^-distributions 
for other models. Even the strong-field approximations ( p^ 
do not agree with each other, and the measure W here im- 
plies a non-constant mean angular velocity in a. 

As we saw above, the overall growth of the deflection 
measure was similar for our multigrid models and for the 
P^M-models with rather large softening parameters (Fig. |^. 
The comparison of the distributions of the velocity deflec- 
tion measure W at the end of the simulations, a = 9 for the 
P^M-models with different softening is shown in Fig. |ll] 
Only the standard model P3M1 (e = 0.2) shows a strong 
collision signature, for the smooth model P3M2 (e = 1) the 
right maximum is already very weak and the 'oversmoothed' 




log W 

Figure 10. Distribution of the velocity defiection measure W 
for the tree-code simulation at different times. The distribution 
is rather wide, but has only one maximum at all times.. 



P3M3 (e = 2) shows only mean-fleld deflections. It is even 
smoother than the APM-model, where the local pairwise in- 
teractions were ignored. We see also that the models with 
a large softening develop much smaller accumulated deflec- 
tions. 
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Figure 11. Distribution of the velocity deflection measure W 
at the flnal stage of simulations for P-^M-models with different 
softening parameters. The role of collisions (the right maximum) 
decreases rapidly with a growing e. 
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Figure 13. Distributions of the deflection rate U for the P'^M- 
model P3M1 for different time intervals. The left maximum de- 
scribes the mean-fleld effects while the right maximum describing 
collisions grows steadily. 




log m 

Figure 12. Distribution of the mean angular velocity lo for dif- 
ferent moments (the model AMGl). The initially concentrated 
distribution spreads in time both due to collisions and mean-field 
effects. 



The deflection characteristics we have studied so far 
have all been accumulated from the start of the simulation. 
As the deflection angle is positively defined, the change of its 
mean value, describing smooth-field effects, will contribute 
in some extent to the velocity deflection parameter W . We 
show a typical distribution of ui in Fig. h2 for the model 



AMGl. Its time dependence is rather slow, but uj tends to 
grow, and this may influence the accumulated velocity de- 
flections. 

In order to eliminate this effect, we performed another 
step of 'differentiation', separating the overall evolution into 
a number of time intervals a £ [a^, Oi+i] and calculating the 
rate of growth of the velocity deflection measure (deflection 
rate) U by 

U^W{ai,ai+Y)/{ai+-i-ai), (13) 

where W{ai,aiJ^\) is the same expression as in ( pT|jl^ ) 
but we have changed the integration limits from (ao,a) to 
{ai,ai+\). It is important to keep in mind that this opera- 
tion does not substract distributions at different times, we 
differentiate along the accumulated angular velocity deflec- 
tions. 

In Fig. ^ we see again the results for the P3M1- 
simulation. The picture is essentially the same as in Fig- 
ure ^ Only the distributions for different times are more 
similar, and the roles of small deflections for a subinterval 
(mean-field effects) and those of large deflections are more 
clearly separated. The small-scale peaks do not move, indi- 
cating that the mean-field effects stay the same during the 
evolution. We built a similar graph for the model P3M4 that 
differs from the present model only by a larger number of 
particles. As suggested by a similar behaviour of the deflec- 
tion measure S (see Fig. ^ the distributions are practically 
the same. 

The distributions of the deflection rate in the adap- 
tive multigrid code are shown in Fig. ^ The figure reveals 
a noticable difference between the [/-distribution from the 
first and the successive time intervals. A probable reason 
for this could be a non-self-consistent initial state that gives 
rise to initial transients (such transients were observed also 
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Figure 14. Distributions of tlie deflection rate U for tlie multi- 
grid model AMGl. After the initial fast evolution towards a 
quasiequilibrium qrowth, most particles keep moving smoothly 
and only a small percentage of them suffers deflections that place 
them in the high-deflection tail of the distribution. 



by Hernquist & Barnes, 1990, in their analysis of tree-code 
models) . This could be caused by the presence of softening 
in the force law as supposed by Hernquist & Barnes. Models 
evolve rapidly away from this stage, but they acquire in the 
process a high-deflection tail of the distribution of the deflec- 
tion measures. The later evolution is much more quiet, with 
most of the deflections coming from small-amplitude mean- 
field accelerations. At later times yet when higher resolution 
grids are used we see also a growth of the high-deflection tail. 
This could be partly a mean-field effect and could be partly 
caused by force errors. This tail is, however, much lower 
than that for the P3Ml-model and there is no sign of any 
maximum. This figure also shows that the strong difference 
between the W-distribution for nonlinear stages (see Fig. ^ 
is due to a large difference between the moments they have 
been constructed, the deflection rate being practically the 
same at a = 9 and at a = 20. 

It is instructive to compare Fig. |l| with Fig. ^, where 
the [/-distribution from the APM simulation is shown. This 
shows the difference between trajectories of particles which 
undergo short-range accelerations and which do not. The 
APM distribution is in fact very close to that for AMG, 
but it almost does not have the long high-defiection tail. 
There are at least two reasons for this - firstly, the adaptive 
grids in the APM did not reach as deep as they did in the 
AMG-models, and secondly, the subgrid force calculation 
procedure used in the AP^M-approach might give cleaner 
forces. 

Extra gravitational collisions inherent to a A^-body code 
can infiuence directly the conclusions that we make com- 
paring our models to observations. As they generate ad- 
ditional velocity changes, their presence is reflected most 
clearly in the velocity distributions. These can be compared 
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Figure 15. Distributions of the deflection rate U for the APM- 
model. Here one observes only the initial fast relaxation, but there 
are practically no collision events. 



easily with observations and they infiuence the properties of 
the gravitationally bound systems that are born during the 
simulations. As an example, we analyzed the distributions 
of relative velocities for particle pairs for different codes. 

The relative velocity vector between a pair of particles 
can be decomposed into the line-of-sight ('radial') compo- 
nent and into the component perpendicular to this direction 
('tangential') (see Gelb 1992). The radial component U|| can 
be defined as 



{vi - V2) ■ {xi - X2) 

'"W = 1 \ 

\Xi — X2\ 

and the perpendicular velocity component v± as 



'v± = {{vi - V2) 



,,2^1/2 



(14) 



(15) 



The dispersions of these velocities are shown in Fig. |l6| 
for the final moment a = 9, corresponding roughly to the 
present epoch, for three different models, P3M1, AMGl and 
APM. This figure shows that the P^M-code gives velocity 
dispersions that are about 50% larger than we get using 
smooth-field codes to model the same patch of the Uni- 
verse. The results we have seen above let us suspect that 
this difference may be mainly due to two-body effects in the 
P'^M-code. And as we already mentioned, this could affect 
the conclusions one makes about the rotation velocities and 
masses of the systems that form and about the ease they 
form with. 



4 CONCLUSIONS 

We have compared a range of P^M-models with different 
softening parameters and spectra and smooth-density mod- 
els of comparable spatial resolution (APM, multigrid). We 
have seen that smooth-field PM-codes are considerably less 
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Figure 16. Dependence of the pairwise longitudinal ("radial") 
velocity dispersion (panel a) and the transversal velocity 

dispersion ("rotation velocities") (7{vj_) on the width of particle 
pairs r for different models. We see that the dispersion obtained 
from the P^M-models is much larger than that from the smooth 
models, even at relatively large distances. 
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colhsional than the P^M-codes. We have found that even 
choosing comoving softening parameters instead of physical 
ones there is a considerable amount of gravitational colli- 
sions in standard P'^M-models. If their influence is crucial 
for the problem the simulation is run for (e.g. the study of 
velocity dispersions or of the formation of bound systems), 
we would recommend to use comoving softening parameters 
e >= 1.0. 

We have also proposed a new approach to study grav- 
itational collisions in nonstationary systems that are com- 
mon in cosmological simulations. We use measures based on 
the accumulated deflection angles of particle orbits and on 
accumulated angular velocity deflections. These measures 
are similar to the commonly used energy diffusion and or- 
bit divergence measures, but they do not saturate during 
evolution. 
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